ExportGridIntegerToNetCDF Subroutine

public subroutine ExportGridIntegerToNetCDF(grid, file, name, action)

Uses

export grid into netcdf file. Two actions are possible:

  • add: add a non-record variable to a netCDF dataset. If file does not exist it is created new. The added variable must have the same dimensions of already present variables. If you try to add a variable already present in netCDF dataset an error is returned.
  • append: append more data to record variables along the unlimited dimension. If netCDF dataset is created new, dimensions and attributes are set, otherwise only new data are written. If current time is already present in netcdf dataset an error is returned

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: grid

grid to be exported
character(len=*), intent(in) :: file

netcdf file to export to
character(len=*), intent(in) :: name

name of variable in netcdf
character(len=*), intent(in) :: action

add or append

Variables

Type Visibility Attributes Name Initial
integer(kind=short), public :: dimTid

id of time dimension
integer(kind=short), public :: dimXid

id of x dimension
integer(kind=short), public :: dimYid

id of y dimension
character(len=7), public :: dummyString
logical, public :: fileExist

true if netcdf exists
integer(kind=short), public :: i

loop index
real(kind=float), public, ALLOCATABLE :: lats(:)

array containing coordinate
real(kind=float), public, ALLOCATABLE :: lons(:)

array containing coordinate
integer(kind=short), public :: mapId

id of grid mapping variable
integer(kind=short), public :: ncId

NetCdf Id for the file
integer(kind=short), public :: ncStatus

error code return by NetCDF routines
integer(kind=short), public :: records

number of records stored in netcdf dataset
type(DateTime), public :: ref_time

reference time to calculate time span
integer(kind=short), public :: slice(3)

used to write new record in right position
character(len=25), public :: string
integer(kind=long), public, ALLOCATABLE :: tempGrid(:,:)
logical, public :: timeExists
integer(kind=short), public :: timeRecord(1)

to write new time record
real(kind=float), public :: timeSpan

time between current time and reference time
real(kind=float), public, ALLOCATABLE :: times(:)

values stored in time variable
integer(kind=short), public, ALLOCATABLE :: varDim(:)

define dimension of variable
logical, public :: varExist

true if variable exists in dataset

integer(kind=short), public :: varId

id of grid variable
integer(kind=short), public :: varTid

id of time variable
integer(kind=short), public :: varXid

id of x variable
integer(kind=short), public :: varYid

id of y variable

Source Code

SUBROUTINE ExportGridIntegerToNetCDF &
!
(grid, file, name, action)

USE StringManipulation, ONLY: &
! imported routines:
StringToUpper, ToString

IMPLICIT NONE

! Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid !! grid to be exported
CHARACTER (LEN = *), INTENT(IN) :: file !!netcdf file to export to
CHARACTER (LEN = *), INTENT(IN) :: name !!name of variable in netcdf
CHARACTER (LEN = *), INTENT(IN) :: action !!add or append

! Local variables:
INTEGER (KIND = short)  :: ncStatus !!error code return by NetCDF routines
INTEGER (KIND = short)  :: ncId  !!NetCdf Id for the file
INTEGER (KIND = short)  :: dimXid  !!id of x dimension
INTEGER (KIND = short)  :: dimYid  !!id of y dimension
INTEGER (KIND = short)  :: dimTid  !!id of time dimension
INTEGER (KIND = short)  :: varXid  !!id of x variable
INTEGER (KIND = short)  :: varYid  !!id of y variable
INTEGER (KIND = short)  :: varTid  !!id of time variable
INTEGER (KIND = short)  :: mapId   !!id of grid mapping variable
INTEGER (KIND = short), ALLOCATABLE  :: varDim (:)  !!define dimension of variable
INTEGER (KIND = short)  :: varId  !!id of grid variable
CHARACTER (LEN = 25)    :: string
REAL (KIND = float), ALLOCATABLE :: lats(:), lons(:) !!array containing coordinate
INTEGER (KIND = short)  :: i !!loop index
LOGICAL :: fileExist    !!true if netcdf exists
LOGICAL :: varExist     !!true if variable exists in dataset 
TYPE (dateTime) :: ref_time !!reference time to calculate time span
CHARACTER (LEN = 7) :: dummyString 
REAL (KIND = float) :: timeSpan !!time between current time and reference time
INTEGER (KIND = short) :: records !!number of records stored in netcdf dataset
INTEGER (KIND = short) :: slice (3) !!used to write new record in right position
INTEGER (KIND = short) :: timeRecord (1) !!to write new time record
REAL (KIND = float), ALLOCATABLE :: times (:) !!values stored in time variable
LOGICAL :: timeExists
INTEGER (KIND = long), ALLOCATABLE :: tempGrid (:,:)

!------------end of declaration------------------------------------------------

!verify if file exists
ncStatus = nf90_open (file, NF90_WRITE, ncId)
IF (ncStatus == nf90_noerr) THEN
  fileExist = .TRUE.
ELSE
  fileExist = .FALSE.
END IF 

!if file does not exist create the file
IF (.NOT. fileExist) THEN
  ncStatus = nf90_create ( file, NF90_CLOBBER, ncId )
  CALL ncErrorHandler (ncStatus)
  
  !define dimensions
  ncStatus = nf90_def_dim ( ncId, 'x', grid % jdim, dimXid )
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_def_dim ( ncId, 'y', grid % idim, dimYid )
  CALL ncErrorHandler (ncStatus)
  !if record (unlimited) variable define time dimension
  IF ( StringToUpper(action) ==  'APPEND' ) THEN
    ncStatus = nf90_def_dim ( ncId, 'time', NF90_UNLIMITED, dimTid )
    CALL ncErrorHandler (ncStatus) 
  END IF
  
  !define coordinate variables
  ncStatus = nf90_def_var (ncId, 'x', NF90_FLOAT, dimXid, varXid)
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_def_var (ncId, 'y', NF90_FLOAT, dimYid, varYid)
  CALL ncErrorHandler (ncStatus)
  !if record (unlimited) variable define time variable
  IF ( StringToUpper(action) ==  'APPEND' ) THEN
    ncStatus = nf90_def_var ( ncId, 'time', NF90_FLOAT, dimTid, varTid )
    CALL ncErrorHandler (ncStatus) 
  END IF
  
  !assign attributes to coordinate variables
  IF (grid % grid_mapping % system == GEODETIC) THEN
      ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'longitude (centre of grid cell)')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'longitude')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'units', 'deg')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'latitude (centre of grid cell)')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'latitude')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'units', 'deg')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y')
      CALL ncErrorHandler (ncStatus)
  ELSE 
      ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'x coordinate of projection')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'projection_x_coordinate')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'units', 'm')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'y coordinate of projection')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'projection_y_coordinate')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'units', 'm')
      CALL ncErrorHandler (ncStatus)
      ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y')
      CALL ncErrorHandler (ncStatus)
END IF

  !if record (unlimited) variable define attributes of time variable
  IF ( StringToUpper(action) ==  'APPEND' ) THEN
    ncStatus = nf90_put_att ( ncId, varTid, 'long_name', 'time' )
    CALL ncErrorHandler (ncStatus)
    ncStatus = nf90_put_att ( ncId, varTid, 'standard_name', 'time' )
    CALL ncErrorHandler (ncStatus)
    string = grid % current_time
    ncStatus = nf90_put_att ( ncId, varTid, 'units', 'seconds since ' // string )
    CALL ncErrorHandler (ncStatus) 
  END IF
  
  !define dummy grid mapping variable
  ncStatus = nf90_def_var (ncId, 'crs', NF90_INT, mapId)
  CALL ncErrorHandler (ncStatus)
  
  !assign attributes to grid mapping variable
  ncStatus = nf90_put_att (ncId, mapId, 'grid_mapping_name', &
                           grid % grid_mapping % name)
  CALL ncErrorHandler (ncStatus)
  DO i = 1, grid % grid_mapping % basic
    ncStatus = nf90_put_att (ncId, mapId, grid % grid_mapping % description (i),&
                            grid % grid_mapping % param (i))
    CALL ncErrorHandler (ncStatus)
  END DO
  !datum
  ncStatus = nf90_put_att (ncId, mapId, 'datum_code', &
                           grid % grid_mapping % datum % code)
  CALL ncErrorHandler (ncStatus)
  
  
  !assign global attributes
  ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'Conventions', 'CF-1.0')
  CALL ncErrorHandler (ncStatus)
  string = UtcNow ()
  ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'history', 'creation date: ' // string )
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'institution', 'Politecnico di Milano' )
  CALL ncErrorHandler (ncStatus)

  !end define mode
  ncStatus = nf90_enddef (ncId)
  CALL ncErrorHandler (ncStatus)

END IF

!enter define mode
ncStatus = nf90_redef (ncId)
CALL ncErrorHandler (ncStatus)
  
!define  variable for the grid
 !If file already exists, have to read dimXd and dimYd
 !and dimTid if necessary
IF ( fileExist) THEN
  ncStatus = nf90_inq_dimid (ncId, "x", dimXid)
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_inq_dimid (ncId, "y", dimYid)
  CALL ncErrorHandler (ncStatus)
  IF (StringToUpper (action) == 'APPEND') THEN
    ncStatus = nf90_inq_dimid (ncId, "time", dimTid)
    CALL ncErrorHandler (ncStatus)
    ncStatus = nf90_inq_varid (ncId, 'time', varTid)
    CALL ncErrorHandler (ncStatus)
  END IF
END IF

!verify if variable already exists 
varId = 0
ncStatus = nf90_inq_varid (ncId, name, varId)
IF (ncStatus == nf90_noerr) THEN
  varExist = .TRUE.
  !if variable already exists in non record dataset: error
  IF (StringToUpper (action) == 'ADD') THEN
    CALL Catch ('error', 'GridLib', 'variable ' // TRIM(name) // &
                ' already exists in netCDf dataset: ', &
                code = ncIOError, argument = file )
  END IF
ELSE
  varExist = .FALSE.

  !define new variable 
  IF (StringToUpper (action) == 'ADD' ) THEN
    ALLOCATE ( varDim (2) )  
  ELSE IF (StringToUpper (action) == 'APPEND' ) THEN
    ALLOCATE ( varDim (3) )
    varDim(3) = dimTid
  END IF
  varDim(1) = dimXid
  varDim(2) = dimYid
  ncStatus = nf90_def_var (ncId, name, NF90_FLOAT, varDim, varId)
  CALL ncErrorHandler (ncStatus)
  DEALLOCATE (varDim)

  !assign attributes to variable
  ncStatus = nf90_put_att (ncId, varId, 'long_name', grid % long_name)
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_put_att (ncId, varId, 'standard_name', grid % standard_name)
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_put_att (ncId, varId, 'units', grid % units)
  CALL ncErrorHandler (ncStatus)
  IF (grid % valid_min /= grid % nodata) THEN
    ncStatus = nf90_put_att (ncId, varId, 'valid_min', grid % valid_min)
    CALL ncErrorHandler (ncStatus)
  END IF
  IF (grid % valid_max /= grid % nodata) THEN
    ncStatus = nf90_put_att (ncId, varId, 'valid_max', grid % valid_max)
    CALL ncErrorHandler (ncStatus)
  END IF
  ncStatus = nf90_put_att (ncId, varId, '_FillValue', grid % nodata)
  CALL ncErrorHandler (ncStatus)
  IF (grid % esri_pe_string /= '') THEN
    ncStatus = nf90_put_att (ncId, varId, 'esri_pe_string', grid % esri_pe_string)
    CALL ncErrorHandler (ncStatus)
  END IF
  ncStatus = nf90_put_att (ncId, varId, 'varying_mode', grid % varying_mode)
  CALL ncErrorHandler (ncStatus)
  ncStatus = nf90_put_att (ncId, varId, 'grid_mapping', 'crs')
  CALL ncErrorHandler (ncStatus)
  
END IF

!end define mode
  ncStatus = nf90_enddef (ncId)
  CALL ncErrorHandler (ncStatus)

IF ( .NOT. fileExist) THEN !if file was created new, add coordinate variables
  !put coordinate
  ALLOCATE ( lons ( grid % jdim) )
  DO i = 1, grid % jdim
    lons (i) = grid % xllcorner + i * grid % cellsize - grid % cellsize / 2.
  END DO
  ncStatus = nf90_put_var (ncId, varXid, lons )
  CALL ncErrorHandler (ncStatus)
  DEALLOCATE ( lons )

  ALLOCATE ( lats ( grid % idim) )
  DO i = 1, grid % idim
    lats (i) = grid % yllcorner + i * grid % cellsize - grid % cellsize / 2.
  END DO
  ncStatus = nf90_put_var (ncId, varYid, lats )
  CALL ncErrorHandler (ncStatus)
  DEALLOCATE ( lats )

END IF

!if unlimited variable put time and grid
IF ( StringToUpper (action) == 'APPEND' ) THEN
  !time
  CALL ParseTime (ncId, ref_time, dummyString)
  timeSpan = grid % current_time - ref_time
  !retrieve number of records already stored in dataset
  ncStatus = nf90_inquire_dimension (ncId, dimTid, len = records)
  CALL ncErrorHandler (ncStatus) 
  !check if time span already exists 
  ALLOCATE ( times (records) )
  ncStatus = nf90_get_var (ncId, varTid, times)
  timeExists = .FALSE.
  DO i = 1, records
    IF ( timeSpan == times (i) ) THEN
      timeExists = .TRUE.
      timeRecord = i
      !(ADD SUPPORT FOR MULTIPLE TIME AXIS
      !CALL Catch ('error', 'GridLib',        &
      !     'time is already present in netCDF dataset: ',  &
      !     code = ncIOError, argument = ToString(timeSpan) )
    END IF
  END DO
  DEALLOCATE (times)
  
  !put time span
  IF ( .NOT. timeExists ) THEN
    timeRecord = records + 1
    ncStatus = nf90_put_var (ncId, varTid, timeSpan, start = timeRecord )
  END IF
  
  !put grid variable
  slice = (/ 1, 1, timeRecord /)
  !swap grid before write to netcdf
  ALLOCATE (tempGrid (grid % jdim, grid % idim) )
  CALL SwapGridIntegerBack (grid % mat, tempGrid) 
  ncStatus = nf90_put_var (ncId, varId, tempGrid, start = slice )
  CALL ncErrorHandler (ncStatus) 
  DEALLOCATE (tempGrid)
ELSE IF (StringToUpper (action) == 'ADD' ) THEN
  !swap grid before write to netcdf
  ALLOCATE (tempGrid (grid % jdim, grid % idim) )
  CALL SwapGridIntegerBack (grid % mat, tempGrid)
  !put variable
  ncStatus = nf90_put_var (ncId, varId, tempGrid )
  CALL ncErrorHandler (ncStatus)
  DEALLOCATE (tempGrid)
END IF

! close file
ncStatus = nf90_close (ncid)
CALL ncErrorHandler (ncStatus)

END SUBROUTINE ExportGridIntegerToNetCDF